Project Description: Examining the enteric bacterial microbiome of West Nile Virus infected mice (pre and post infection) to determine alterations in bacterial community strucutre or specific taxa which covary with lethality.
Primary Collaborator: Larissa Thakray (lthackray@wustl.edu)
Other relevant files: The following Phyloseq objects are available. Each is distinguished based on the 16S reference database used for taxonomic classification. RDP and Silva were processed through the species assignment workflow:
References: * http://f1000research.com/articles/5-1492/v1
Workflow details: The R commands below represent a full analysis of the following:
library("tidyverse")
packageVersion("tidyverse")
## [1] '1.2.1'
library("reshape2")
packageVersion("reshape2")
## [1] '1.4.3'
library("plyr")
packageVersion("plyr")
## [1] '1.8.4'
library("phyloseq")
packageVersion("phyloseq")
## [1] '1.23.1'
library("RColorBrewer")
packageVersion("RColorBrewer")
## [1] '1.1.2'
library("vegan")
packageVersion("vegan")
## [1] '2.4.5'
library("gridExtra")
packageVersion("gridExtra")
## [1] '2.3'
library("knitr")
packageVersion("knitr")
## [1] '1.18'
library("plotly")
packageVersion("plotly")
## [1] '4.7.1'
library("microbiome")
packageVersion("microbiome")
## [1] '1.0.0'
library("ggpubr")
packageVersion("ggpubr")
## [1] '0.1.6'
library("data.table")
packageVersion("data.table")
## [1] '1.10.4.3'
library("mgcv")
packageVersion("mgcv")
## [1] '1.8.22'
library("pairwiseAdonis") #https://github.com/pmartinezarbizu/pairwiseAdonis
packageVersion("pairwiseAdonis")
## [1] '0.0.1'
# Load Phyloseq Object
# Selected RDP due to it's up-to-date nature and conservative taxonomy. Other files are also valid for anlysis but are not explored here
ps0 <- readRDS("./Data/ps0.wnv_antibiotics.rdp.RDS")
ps0
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3229 taxa and 535 samples ]
## tax_table() Taxonomy Table: [ 3229 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3229 tips and 3227 internal nodes ]
# Load mapping file
map <- import_qiime_sample_data("./Data/mapping_wnv_antibiotics.txt")
dim(map)
## [1] 520 24
ps0 <- merge_phyloseq(ps0, map)
ps0
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3229 taxa and 520 samples ]
## sample_data() Sample Data: [ 520 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 3229 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3229 tips and 3227 internal nodes ]
# Perform a few sanity checks
sample_variables(ps0)
## [1] "X.SampleID" "BarcodeSequence" "LinkerPrimerSequence"
## [4] "Name" "MiSeq" "DNAplate"
## [7] "Well" "Sample.ID" "Description"
## [10] "X16SPlate" "SampleNo" "Run"
## [13] "DNAcell" "Sex" "CageBeforeTreatment"
## [16] "MouseBeforeTreatment" "CageDuringTreatment" "DaysTreatment"
## [19] "Day" "Treatment" "Virus"
## [22] "SurvivalStatus" "DayDeadPostInfection" "Description.1"
ntaxa(ps0)
## [1] 3229
rank_names(ps0)
## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"
get_taxa_unique(ps0, "Phylum")
## [1] "Bacteroidetes" "Proteobacteria"
## [3] "Firmicutes" "Verrucomicrobia"
## [5] "Tenericutes" NA
## [7] "Actinobacteria" "Cyanobacteria/Chloroplast"
## [9] "Fusobacteria"
# Remove Day -14 cohoused data
levels(sample_data(ps0)$DaysTreatment)
## [1] "D.14" "D0" "D13" "D16" "D18" "D20" "D3" "D7"
ps0 <- subset_samples(ps0, DaysTreatment != "D.14")
levels(sample_data(ps0)$DaysTreatment)
## [1] "D0" "D13" "D16" "D18" "D20" "D3" "D7"
# Remove uninfected samples
levels(sample_data(ps0)$Virus)
## [1] "Uninfected" "WNV2000"
ps0 <- subset_samples(ps0, Virus == "WNV2000")
levels(sample_data(ps0)$Virus)
## [1] "WNV2000"
# Remove taxa no longer part of the count table due to sample removal
summary(taxa_sums(ps0))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 3 2353 10 928695
ps0 <- prune_taxa(taxa_sums(ps0) > 0, ps0)
summary(taxa_sums(ps0))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.0 4.0 9.0 4233.2 72.5 928695.0
# Factor re-ordering, relabelling, etc.
# Reorder Time points
levels(sample_data(ps0)$DaysTreatment)
## [1] "D0" "D13" "D16" "D18" "D20" "D3" "D7"
sample_data(ps0)$DaysTreatment <- factor(sample_data(ps0)$DaysTreatment, levels = c("D0", "D3", "D7", "D13", "D16", "D18", "D20"))
levels(sample_data(ps0)$DaysTreatment)
## [1] "D0" "D3" "D7" "D13" "D16" "D18" "D20"
# Reorder Treatments
levels(sample_data(ps0)$Treatment)
## [1] "Amp" "AmpMetro" "Metro" "Vehicle"
sample_data(ps0)$Treatment <- factor(sample_data(ps0)$Treatment, levels = c("Vehicle","Metro","Amp","AmpMetro"))
levels(sample_data(ps0)$Treatment)
## [1] "Vehicle" "Metro" "Amp" "AmpMetro"
# Relabel Treatments
sample_data(ps0)$Treatment <- factor(sample_data(ps0)$Treatment, labels = c("Vehicle","Metro","Amp","Amp + Metro"))
levels(sample_data(ps0)$Treatment)
## [1] "Vehicle" "Metro" "Amp" "Amp + Metro"
# Create a new data frame of the sorted row sums, a column of sorted values from 1 to the total number of individuals/counts for each ASV and a categorical variable stating these are all ASVs.
readsumsdf = data.frame(nreads = sort(taxa_sums(ps0), TRUE),
sorted = 1:ntaxa(ps0),
type = "ASVs")
# Add a column of sample sums (total number of individuals per sample)
readsumsdf = rbind(readsumsdf,
data.frame(nreads = sort(sample_sums(ps0), TRUE),
sorted = 1:nsamples(ps0),
type = "Samples"))
# Make a data frame with a column for the read counts of each sample for histogram production
sample_sum_df <- data.frame(sum = sample_sums(ps0))
# Make plots
# Generates a bar plot with # of reads (y-axis) for each taxa. Sorted from most to least abundant
# Generates a second bar plot with # of reads (y-axis) per sample. Sorted from most to least
p.reads = ggplot(readsumsdf, aes(x = sorted, y = nreads)) +
geom_bar(stat = "identity") +
ggtitle("ASV Assessment") +
scale_y_log10() +
facet_wrap(~type, scales = "free") +
ylab("# of Reads")
# Histogram of the number of Samples (y-axis) at various read depths
p.reads.hist <- ggplot(sample_sum_df, aes(x = sum)) +
geom_histogram(color = "black", fill = "firebrick3", binwidth = 150) +
ggtitle("Distribution of sample sequencing depth") +
xlab("Read counts") +
ylab("# of Samples")
# Final plot, side-by-side
grid.arrange(p.reads, p.reads.hist, ncol = 2)
# Basic summary statistics
summary(sample_sums(ps0))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 86 2064 19679 24122 37592 97513
# Format a data table to combine sample summary data with sample variable data
ss <- sample_sums(ps0)
sd <- as.data.frame(sample_data(ps0))
ss.df <- merge(sd, data.frame("ASVs" = ss), by ="row.names")
# Plot the data by the treatment variable
y = 1000 # Set a threshold for the minimum number of acceptable reads. Can start as a guess
x = "DaysTreatment" # Set the x-axis variable you want to examine
label = "Sample.ID" # This is the label you want to overlay on the points that are below threshold y. Should be something sample specific
p.ss.boxplot <- ggplot(ss.df, aes_string(x, y = "ASVs")) +
stat_boxplot(geom = "errorbar", position = position_dodge(width = 0.8)) +
geom_boxplot(outlier.colour="NA", position = position_dodge(width = 0.8), alpha = 0.2) +
scale_y_log10() +
facet_wrap(~Treatment) +
geom_hline(yintercept = y, lty = 2) +
geom_point(position=position_jitterdodge(dodge.width = 0.8), aes_string(color = "SurvivalStatus"), size = 1.2) +
geom_text(data = ss.df, aes_string(x, y="ASVs", label=label), size=2) # This labels a subset that fall below threshold variable y and labels them with the label variable
p.ss.boxplot
write.table(ss.df, file = "./Results/asv_stats.txt", sep = "\t")
List of samples selected as outliers: c(“D16.M5”, “D16.M2”, “D18.M3”, “D18.M5”, “D18.K1”, “D20.M5”, “D20.M2”)
# Outlier samples: c("D16.M5", "D16.M2", "D18.M3", "D18.M5", "D18.K1", "D20.M5", "D20.M2")
nsamples(ps0)
## [1] 315
ps1 <- ps0 %>%
subset_samples(
Sample.ID != "D16.M5" &
Sample.ID != "D16.M2" &
Sample.ID != "D18.M3" &
Sample.ID != "D18.M5" &
Sample.ID != "D18.K1" &
Sample.ID != "D20.M5" &
Sample.ID != "D20.M2"
)
nsamples(ps1)
## [1] 309
saveRDS(ps1, file = "./Results/ps1.RDS")
ggpaired(ss.df, x = "DaysTreatment", y = "ASVs", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "gray") +
scale_y_log10() +
facet_grid(~Treatment) +
stat_compare_means(ref.group = "D0", hide.ns = TRUE, label = "p.signif") +
theme(axis.text.x = element_text(size = 8)) +
theme(axis.text.y = element_text(size = 8)) +
theme(axis.title.x = element_blank()) +
ylab("Read Counts")
##Overall sample relationship to evaluate sample outliers
# Begin by removing sequences that were classified as either mitochondria or chlorplast
ntaxa(ps1) # Check the number of taxa prior to removal
## [1] 1795
ps1 <- ps1 %>%
subset_taxa(
Family != "mitochondria" &
Class != "Chloroplast"
)
ntaxa(ps1) # Confirm that the taxa were removed
## [1] 1372
# Create a data table for ggploting
ps1_phylum <- ps1 %>%
tax_glom(taxrank = "Phylum") %>% # agglomerate at phylum level
transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance (or use ps0.ra)
psmelt() %>% # Melt to long format for easy ggploting
filter(Abundance > 0.01) # Filter out low abundance taxa
# Convert Sample No to a factor because R is weird sometime
ps1_phylum$SampleNo <- as.factor(ps1_phylum$SampleNo)
# Plot - Phylum
p.ra.phylum <- ggplot(ps1_phylum, aes(x = SampleNo, y = Abundance, fill = Phylum)) +
geom_bar(stat = "identity", width = 1) +
facet_wrap(Treatment~DaysTreatment, scales = "free_x", nrow = 4, ncol = 7) +
theme(axis.text.x = element_blank()) +
theme(axis.title.x = element_blank()) +
ggtitle("Abundant Phylum (> 1%)")
p.ra.phylum
# Note: This is a nice place to output tables of data that you may want to use for other analysis, or to include as supplemental data for publication
# You can rerun the first bit of code in this chunk and change Phylum to Species for a table with all possible classifications
write.table(ps1_phylum, file = "./Results/phylum_relab.txt", sep = "\t")
ggplotly(p.ra.phylum)
# agglomerate taxa
glom <- tax_glom(ps1.ra, taxrank = 'Phylum')
# create dataframe from phyloseq object
dat <- as.tibble(psmelt(glom))
# convert Phylum to a character vector from a factor because R
# dat$Phylum <- as.character(dat$Phylum)
# Reorder Phylum levels from most -> least abundant
levels(dat$Phylum)
## [1] "Actinobacteria" "Bacteroidetes" "Firmicutes" "Fusobacteria"
## [5] "Proteobacteria" "Tenericutes" "Verrucomicrobia"
dat$Phylum <- factor(dat$Phylum, levels = c("Bacteroidetes", "Firmicutes", "Proteobacteria", "Tenericutes", "Actinobacteria", "Verrucomicrobia"))
levels(dat$Phylum)
## [1] "Bacteroidetes" "Firmicutes" "Proteobacteria" "Tenericutes"
## [5] "Actinobacteria" "Verrucomicrobia"
levels(dat$Treatment)
## [1] "Vehicle" "Metro" "Amp" "Amp + Metro"
dat$Treatment <- factor(dat$Treatment, levels = c("Vehicle", "Metro", "Amp", "Amp + Metro"))
levels(dat$Treatment)
## [1] "Vehicle" "Metro" "Amp" "Amp + Metro"
p.boxplot.phylum.1 <- ggpaired(dat, x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
#geom_point(size = 1) +
#geom_jitter(width = 0.2, alpha = 0.7) +
ylab("Relative Abundance") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
facet_grid(Treatment~Phylum) +
theme(axis.text.x = element_text(size = 6)) +
theme(axis.text.y = element_text(size = 6)) +
theme(axis.title.y = element_text(size = 6)) +
theme(axis.title.x = element_blank()) +
theme(strip.text = element_text(size = 6)) +
ylim(0,1.2)
p.boxplot.phylum.1
# Reduced to most abundant phylum
summary(dat$Phylum)
## Bacteroidetes Firmicutes Proteobacteria Tenericutes
## 309 309 309 309
## Actinobacteria Verrucomicrobia NA's
## 309 309 309
dat.1 <- filter(dat, Phylum %in% c("Bacteroidetes", "Firmicutes", "Proteobacteria", "Tenericutes"))
dat.1 <- droplevels(dat.1)
summary(dat.1$Phylum)
## Bacteroidetes Firmicutes Proteobacteria Tenericutes
## 309 309 309 309
levels(dat.1$Treatment)
## [1] "Vehicle" "Metro" "Amp" "Amp + Metro"
dat.1$Treatment <- factor(dat.1$Treatment, levels = c("Vehicle", "Amp", "Metro", "Amp + Metro"))
levels(dat.1$Treatment)
## [1] "Vehicle" "Amp" "Metro" "Amp + Metro"
p.boxplot.phylum.2 <- ggpaired(dat.1, x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
#geom_point(size = 1) +
#geom_jitter(width = 0.2, alpha = 0.7) +
ylab("Relative Abundance") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 0.95) +
facet_grid(Treatment~Phylum) +
theme(axis.text.x = element_text(size = 6)) +
theme(axis.text.y = element_text(size = 6)) +
theme(axis.title.y = element_text(size = 6)) +
theme(axis.title.x = element_blank()) +
theme(strip.text = element_text(size = 6)) +
ylim(0,1)
p.boxplot.phylum.2
## Phylum level general additive model (GAM) plots
# Define color scheme
my.cols <- brewer.pal(n = 8, "Dark2")
my.cols[3] <- "#08519C"
# Phyla plots with GAM smoother
p.gam.phylum <- ggplot(dat.1, aes(x = Day, y = Abundance, color = Phylum, group = Phylum)) +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7)) +
facet_grid(~Treatment) +
ylab("Relative Abundance") +
geom_point(size = 1.25, alpha = 0.4) +
theme(axis.text.x = element_text(size = 10)) +
theme(axis.text.y = element_text(size = 10)) +
theme(axis.title.y = element_text(size = 10)) +
theme(axis.title.x = element_blank()) +
theme(strip.text = element_text(size = 10)) +
scale_y_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2)) +
scale_x_continuous(breaks = c(0, 3, 7, 13, 16, 18, 20)) +
scale_color_manual(values = my.cols) +
theme(strip.background = element_blank()) +
theme(strip.text.x = element_blank()) +
theme(axis.title.y = element_blank()) +
guides(color=guide_legend(override.aes=list(fill=NA))) +
theme(legend.title = element_blank())
p.gam.phylum
# Testing plots - Close up view of how each phyla behaves
dat.bac <- subset(dat, Phylum == "Bacteroidetes")
dat.firm <- subset(dat, Phylum == "Firmicutes")
dat.teneri <- subset(dat, Phylum == "Tenericutes")
dat.proteo <- subset(dat, Phylum == "Proteobacteria")
# Convert data frame so that each animal (CageDuringTreatment) has a unique phyla and abundnace (long form transformation with melt)
dat.1.melt <- melt(dat.1, id.vars = c("CageDuringTreatment", "Day", "Phylum", "Treatment"), measure.vars = c("Abundance"))
# Recast the data by indiviudal phyla
dat.1.cast <- dcast(dat.1.melt, CageDuringTreatment + Treatment + Day ~ Phylum)
# Subset melt/casted frames by Treatment for within treatment testing
dat.1.cast.vehicle <- filter(dat.1.cast, Treatment %in% "Vehicle")
dat.1.cast.amp <- filter(dat.1.cast, Treatment %in% "Amp")
dat.1.cast.metro <- filter(dat.1.cast, Treatment %in% "Metro")
dat.1.cast.ampmetro <- filter(dat.1.cast, Treatment %in% "Amp + Metro")
## General Additive Model (GAM) testing for the contribution of the relative abundance of each Phyla to the change in abundnace over time
# Vehicle
mod_gam.vehicle.1 <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.vehicle)
summary(mod_gam.vehicle.1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr",
## k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes,
## bs = "cr", k = 7)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.000 0.742 14.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Bacteroidetes) 0.9999 1.000 3.823 0.0550 .
## s(Firmicutes) 1.0000 1.000 4.377 0.0404 *
## s(Proteobacteria) 2.0598 2.473 2.946 0.0457 *
## s(Tenericutes) 2.5961 3.102 1.925 0.1226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.261 Deviance explained = 33.3%
## GCV = 43.271 Scale est. = 38.539 n = 70
ggpaired(subset(dat.proteo, Treatment == "Vehicle"), x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
ggtitle("Proteobacteria in Vehicle Control Mice")
# Metro
mod_gam.metro.1 <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.metro)
summary(mod_gam.metro.1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr",
## k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes,
## bs = "cr", k = 7)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.0000 0.6067 18.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Bacteroidetes) 5.467 5.796 1.710 0.1949
## s(Firmicutes) 4.314 5.084 1.582 0.1809
## s(Proteobacteria) 4.694 5.250 2.766 0.0205 *
## s(Tenericutes) 5.809 5.969 8.241 5.4e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.506 Deviance explained = 65.1%
## GCV = 37.023 Scale est. = 25.766 n = 70
ggpaired(subset(dat.teneri, Treatment == "Metro"), x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
ggtitle("Tenericutes in Metro Treated Mice")
ggpaired(subset(dat.proteo, Treatment == "Metro"), x = "DaysTreatment", y = "Abundance", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "lightgray", point.size = 0.3, width = 0.4) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE, label.y.npc = 1) +
ggtitle("Proteobacteria in Metro Treated Mice")
# Amp
mod_gam.amp.1 <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.amp)
summary(mod_gam.amp.1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr",
## k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes,
## bs = "cr", k = 7)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.0000 0.5499 20 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Bacteroidetes) 1.000 1.000 0.193 0.661842
## s(Firmicutes) 5.023 5.596 3.057 0.016504 *
## s(Proteobacteria) 5.264 5.632 3.437 0.009471 **
## s(Tenericutes) 4.725 5.007 4.598 0.000863 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.388 Deviance explained = 48.3%
## GCV = 37.89 Scale est. = 31.751 n = 105
# Amp + Metro
mod_gam.ampmetro <- gam(Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr", k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes, bs = "cr", k = 7), data = dat.1.cast.ampmetro)
summary(mod_gam.ampmetro)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Day ~ s(Bacteroidetes, bs = "cr", k = 7) + s(Firmicutes, bs = "cr",
## k = 7) + s(Proteobacteria, bs = "cr", k = 7) + s(Tenericutes,
## bs = "cr", k = 7)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.3438 0.5071 20.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Bacteroidetes) 1.0544 1.2899 28.445 0.000270 ***
## s(Firmicutes) 0.7196 0.7678 0.065 0.824564
## s(Proteobacteria) 0.4675 0.4717 140.607 3.39e-12 ***
## s(Tenericutes) 2.3332 2.6542 7.372 0.000382 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 23/25
## R-sq.(adj) = 0.682 Deviance explained = 70.5%
## GCV = 18.031 Scale est. = 16.46 n = 64
# Diversity
diversity <- global(ps1)
head(diversity)
## richness_0 richness_20 richness_50 richness_80
## 101.Thackray.D0.E1 121 121 121 121
## 102.Thackray.D0.E2 166 166 166 166
## 103.Thackray.D0.E3 93 93 93 93
## 104.Thackray.D0.E4 151 151 151 151
## 105.Thackray.D0.E5 173 173 173 173
## 106.Thackray.D0.F1 102 102 102 102
## diversities_inverse_simpson diversities_gini_simpson
## 101.Thackray.D0.E1 12.170885 0.9178367
## 102.Thackray.D0.E2 25.401667 0.9606325
## 103.Thackray.D0.E3 10.355095 0.9034292
## 104.Thackray.D0.E4 24.015585 0.9583604
## 105.Thackray.D0.E5 30.875326 0.9676117
## 106.Thackray.D0.F1 9.038626 0.8893637
## diversities_shannon diversities_fisher
## 101.Thackray.D0.E1 3.081963 16.60700
## 102.Thackray.D0.E2 3.983350 23.16563
## 103.Thackray.D0.E3 2.923112 12.47155
## 104.Thackray.D0.E4 3.802921 21.14357
## 105.Thackray.D0.E5 4.068889 24.48143
## 106.Thackray.D0.F1 2.836865 13.94334
## diversities_coverage evenness_camargo evenness_pielou
## 101.Thackray.D0.E1 5 0.01410119 0.6426391
## 102.Thackray.D0.E2 10 0.03354858 0.7792174
## 103.Thackray.D0.E3 4 0.01148292 0.6449084
## 104.Thackray.D0.E4 8 0.02752578 0.7579648
## 105.Thackray.D0.E5 11 0.03541365 0.7895708
## 106.Thackray.D0.F1 3 0.01154458 0.6133797
## evenness_simpson evenness_evar evenness_bulla
## 101.Thackray.D0.E1 0.008929483 0.2255497 0.06979504
## 102.Thackray.D0.E2 0.018636586 0.2804325 0.10587282
## 103.Thackray.D0.E3 0.007597282 0.2026502 0.05605963
## 104.Thackray.D0.E4 0.017619651 0.2691205 0.09552794
## 105.Thackray.D0.E5 0.022652477 0.2584945 0.10681799
## 106.Thackray.D0.F1 0.006631420 0.2367310 0.06282215
## dominance_dbp dominance_dmn dominance_absolute
## 101.Thackray.D0.E1 0.16531143 0.2888926 4005
## 102.Thackray.D0.E2 0.13306855 0.2103331 3987
## 103.Thackray.D0.E3 0.20086163 0.3593830 4336
## 104.Thackray.D0.E4 0.09593916 0.1840488 2561
## 105.Thackray.D0.E5 0.10275191 0.1638590 2946
## 106.Thackray.D0.F1 0.19345906 0.3746479 4052
## dominance_relative dominance_simpson
## 101.Thackray.D0.E1 0.16531143 0.08216330
## 102.Thackray.D0.E2 0.13306855 0.03936749
## 103.Thackray.D0.E3 0.20086163 0.09657082
## 104.Thackray.D0.E4 0.09593916 0.04163963
## 105.Thackray.D0.E5 0.10275191 0.03238832
## 106.Thackray.D0.F1 0.19345906 0.11063629
## dominance_core_abundance dominance_gini
## 101.Thackray.D0.E1 0.7752920 0.9858988
## 102.Thackray.D0.E2 0.5431547 0.9664514
## 103.Thackray.D0.E3 0.7725483 0.9885171
## 104.Thackray.D0.E4 0.6345995 0.9724742
## 105.Thackray.D0.E5 0.4216804 0.9645863
## 106.Thackray.D0.F1 0.8083075 0.9884554
## rarity_log_modulo_skewness rarity_low_abundance
## 101.Thackray.D0.E1 2.059941 0.05733273
## 102.Thackray.D0.E2 2.058710 0.07970095
## 103.Thackray.D0.E3 2.059384 0.04660212
## 104.Thackray.D0.E4 2.059367 0.06799281
## 105.Thackray.D0.E5 2.060132 0.06867567
## 106.Thackray.D0.F1 2.060471 0.05743614
## rarity_noncore_abundance rarity_rare_abundance
## 101.Thackray.D0.E1 0.08527676 0.03178272
## 102.Thackray.D0.E2 0.27635004 0.09892531
## 103.Thackray.D0.E3 0.10779636 0.04377635
## 104.Thackray.D0.E4 0.19850903 0.05952649
## 105.Thackray.D0.E5 0.28607304 0.11677305
## 106.Thackray.D0.F1 0.10188589 0.06378611
ps1.rich <- merge(sd, diversity, by ="row.names") # merge sd.1 by row names
# Add divergence measurements
ps1.rich$divergence <- divergence(ps1)
p.rich.treatment <- ggpaired(ps1.rich, x = "DaysTreatment", y = "richness_0", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "gray") +
geom_jitter(width = 0.2) +
facet_grid(~Treatment) +
ylab("Richness") +
theme(axis.title.x = element_blank()) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE)
p.sd.treatment <- ggpaired(ps1.rich, x = "DaysTreatment", y = "diversities_shannon", outlier.shape = NA, id = "CageDuringTreatment", line.size = 0.5, line.color = "gray") +
geom_boxplot(position = position_dodge()) +
geom_jitter(width = 0.2) +
facet_grid(~Treatment) +
ylab("Shannon diversity") +
theme(axis.title.x = element_blank()) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "D0", hide.ns = TRUE)
grid.arrange(p.rich.treatment, p.sd.treatment, nrow = 2)
## Alpha diversity general additive model (GAM) analysis.
ps1.rich.melt <- melt(ps1.rich, id.vars = c("Treatment", "Day", "DaysTreatment"), measure.vars = c("richness_0"))
ps1.sd.melt <- melt(ps1.rich, id.vars = c("Treatment", "Day", "DaysTreatment"), measure.vars = c("diversities_shannon"))
# Richness
p.rich.gam.treat <- ggplot(ps1.rich.melt, aes(x = Day, y = value, color = Treatment, group = Treatment)) +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7)) +
ylab("Richness") +
geom_point(size = 1.25, alpha = 0.5) +
theme(axis.text.x = element_text(size = 10)) +
theme(axis.text.y = element_text(size = 10)) +
theme(axis.title.y = element_blank()) +
theme(axis.title.x = element_blank()) +
theme(strip.text = element_text(size = 10)) +
scale_color_manual(values = c("black", "chocolate", "green", "purple")) +
theme(legend.position = "NULL") +
scale_y_continuous(limits = c(0,250), breaks = c(0,50, 100, 150, 200, 250)) +
scale_x_continuous(breaks = c(0,3,7,13,16,18,20))
# Shannon diversity
p.sd.gam.treat <- ggplot(ps1.sd.melt, aes(x = Day, y = value, color = Treatment, group = Treatment)) +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7)) +
ylab("Shannon Diversity") +
geom_point(size = 1.25, alpha = 0.5) +
theme(axis.text.x = element_text(size = 10)) +
theme(axis.text.y = element_text(size = 10)) +
theme(axis.title.y = element_blank()) +
theme(axis.title.x = element_blank()) +
theme(strip.text = element_text(size = 10)) +
scale_color_manual(values = c("black", "green", "chocolate", "purple"), labels = c("Vehicle", "A", "M", "AM")) +
guides(color=guide_legend(override.aes=list(fill=NA))) +
theme(legend.position = "right") +
theme(legend.title = element_blank()) +
scale_y_continuous(limits = c(0,5), breaks = c(0,1,2,3,4,5)) +
scale_x_continuous(breaks = c(0,3,7,13,16,18,20))
#theme(legend.title = element_blank())
grid.arrange(p.rich.gam.treat, p.sd.gam.treat, ncol = 2)
## Ordination
#Ordination Analysis
#Beta Diversity has same trend of timepoints with longtail and bimodal read counts having larger elipses
ord.pcoa.bray <- ordinate(ps1, method = "PCoA", distance = "bray")
ord.pcoa.uni <- ordinate(ps1, method = "PCoA", distance = "unifrac")
ord.pcoa.wuni <- ordinate(ps1, method = "PCoA", distance = "wunifrac")
## Ordination plots all samples
# Bray
p.pcoa.bray <- plot_ordination(ps1, ord.pcoa.bray, color = "SurvivalStatus", axes = c(1,2)) +
geom_point(size = 2) +
ggtitle("PCoA of UniFrac Distances") +
facet_grid(Treatment~DaysTreatment)
#stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = SurvivalStatus))
p.pcoa.bray
# Unifrac
p.pcoa.uni <- plot_ordination(ps1, ord.pcoa.uni, color = "SurvivalStatus", axes = c(1,2)) +
geom_point(size = 2) +
ggtitle("PCoA of UniFrac Distances") +
facet_grid(Treatment~DaysTreatment)
#stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = SurvivalStatus))
p.pcoa.uni
# Weighted Unifrac
p.pcoa.wuni <- plot_ordination(ps1, ord.pcoa.wuni, color = "SurvivalStatus", axes = c(1,2)) +
geom_point(size = 2) +
ggtitle("PCoA of wUniFrac Distances") +
facet_grid(Treatment~DaysTreatment) +
stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = SurvivalStatus))
p.pcoa.wuni
p.pcoa.uni.treat <- plot_ordination(ps1, ord.pcoa.uni, color = "Treatment", shape = "SurvivalStatus") +
geom_point(size = 3) +
# ggtitle("PCoA of UniFrac Distances") +
facet_grid(~DaysTreatment) +
scale_color_manual(values = c("black", "green", "chocolate", "purple"), labels = c("Vehicle", "A", "M", "AM")) +
theme(axis.title.y = element_blank()) +
theme(axis.title.x = element_blank()) +
theme(legend.position = "bottom") +
theme(legend.text = element_text(size = 10)) +
theme(legend.title = element_blank()) +
theme(strip.background = element_blank()) +
theme(strip.text.x = element_blank())
p.pcoa.uni.treat
##Group significance testing with ADONIS
# Set a random seed so that exact results can be reproduced
set.seed(10000)
# Function to run adonis test on a physeq object and a variable from metadata
doadonis <- function(physeq, category) {
bdist <- phyloseq::distance(physeq, "unifrac")
col <- as(sample_data(physeq), "data.frame")[ ,category]
# Adonis test
adonis.bdist <- adonis(bdist ~ col)
print("Adonis results:")
print(adonis.bdist)
# Homogeneity of dispersion test
betatax = betadisper(bdist,col)
p = permutest(betatax)
print("Betadisper results:")
print(p$tab)
}
doadonis(ps1, "Treatment")
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 3 17.812 5.9375 39.413 0.27937 0.001 ***
## Residuals 305 45.947 0.1506 0.72063
## Total 308 63.760 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 4.419388 1.473129306 262.2619 999 0.001
## Residuals 305 1.713190 0.005617016 NA NA NA
library(pairwiseAdonis)
ps1.d0 <- subset_samples(ps1, DaysTreatment == "D0")
ps1_otu_table.d0 <- as.data.frame(otu_table(ps1.d0))
sd.df.d0 <- as.data.frame(sample_data(ps1.d0))
ps1_otu_table.d0$Treatment <- sd.df.d0$Treatment
ps1.d3 <- subset_samples(ps1, DaysTreatment == "D3")
ps1_otu_table.d3 <- as.data.frame(otu_table(ps1.d3))
sd.df.d3 <- as.data.frame(sample_data(ps1.d3))
ps1_otu_table.d3$Treatment <- sd.df.d3$Treatment
ps1.d7 <- subset_samples(ps1, DaysTreatment == "D7")
ps1_otu_table.d7 <- as.data.frame(otu_table(ps1.d7))
sd.df.d7 <- as.data.frame(sample_data(ps1.d7))
ps1_otu_table.d7$Treatment <- sd.df.d7$Treatment
ps1.d13 <- subset_samples(ps1, DaysTreatment == "D13")
ps1_otu_table.d13 <- as.data.frame(otu_table(ps1.d13))
sd.df.d13 <- as.data.frame(sample_data(ps1.d13))
ps1_otu_table.d13$Treatment <- sd.df.d13$Treatment
ps1.d16 <- subset_samples(ps1, DaysTreatment == "D16")
ps1_otu_table.d16 <- as.data.frame(otu_table(ps1.d16))
sd.df.d16 <- as.data.frame(sample_data(ps1.d16))
ps1_otu_table.d16$Treatment <- sd.df.d16$Treatment
ps1.d20 <- subset_samples(ps1, DaysTreatment == "D20")
ps1_otu_table.d20 <- as.data.frame(otu_table(ps1.d20))
sd.df.d20 <- as.data.frame(sample_data(ps1.d20))
ps1_otu_table.d20$Treatment <- sd.df.d20$Treatment
pairwise.adonis(ps1_otu_table.d0[,1:1363], ps1_otu_table.d0$Treatment)
## [1] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
## pairs F.Model R2 p.value p.adjusted sig
## 1 Metro vs Amp 2.541514 0.09950522 0.028 0.168
## 2 Metro vs Amp + Metro 3.012322 0.14335980 0.008 0.048 .
## 3 Metro vs Vehicle 3.008221 0.14319254 0.020 0.120
## 4 Amp vs Amp + Metro 0.952873 0.03978116 0.465 1.000
## 5 Amp vs Vehicle 2.675326 0.10419834 0.010 0.060
## 6 Amp + Metro vs Vehicle 3.143583 0.14867786 0.001 0.006 *
pairwise.adonis(ps1_otu_table.d3[,1:1363], ps1_otu_table.d3$Treatment)
## [1] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
## pairs F.Model R2 p.value p.adjusted sig
## 1 Vehicle vs Metro 8.531852 0.3215702 0.001 0.006 *
## 2 Vehicle vs Amp 17.429204 0.4311043 0.001 0.006 *
## 3 Vehicle vs Amp + Metro 8.652032 0.3246294 0.001 0.006 *
## 4 Metro vs Amp 20.688576 0.4735466 0.001 0.006 *
## 5 Metro vs Amp + Metro 11.193323 0.3834207 0.001 0.006 *
## 6 Amp vs Amp + Metro 3.885616 0.1445240 0.001 0.006 *
pairwise.adonis(ps1_otu_table.d7[,1:1363], ps1_otu_table.d7$Treatment)
## [1] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
## pairs F.Model R2 p.value p.adjusted sig
## 1 Vehicle vs Metro 5.590241 0.2369726 0.001 0.006 *
## 2 Vehicle vs Amp 18.023438 0.4393449 0.001 0.006 *
## 3 Vehicle vs Amp + Metro 11.673062 0.3933892 0.001 0.006 *
## 4 Metro vs Amp 27.042017 0.5403862 0.001 0.006 *
## 5 Metro vs Amp + Metro 16.419568 0.4770417 0.001 0.006 *
## 6 Amp vs Amp + Metro 4.396955 0.1604906 0.001 0.006 *
pairwise.adonis(ps1_otu_table.d13[,1:1363], ps1_otu_table.d13$Treatment)
## [1] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
## pairs F.Model R2 p.value p.adjusted sig
## 1 Vehicle vs Metro 3.723724 0.1714128 0.002 0.012 .
## 2 Vehicle vs Amp 19.990937 0.4650035 0.001 0.006 *
## 3 Vehicle vs Amp + Metro 18.657895 0.5089734 0.001 0.006 *
## 4 Metro vs Amp 20.031006 0.4655017 0.001 0.006 *
## 5 Metro vs Amp + Metro 19.103878 0.5148755 0.001 0.006 *
## 6 Amp vs Amp + Metro 9.760540 0.2979359 0.001 0.006 *
pairwise.adonis(ps1_otu_table.d16[,1:1363], ps1_otu_table.d16$Treatment)
## [1] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
## pairs F.Model R2 p.value p.adjusted sig
## 1 Vehicle vs Metro 2.568894 0.1248922 0.014 0.084
## 2 Vehicle vs Amp 20.190996 0.4674816 0.001 0.006 *
## 3 Vehicle vs Amp + Metro 32.948284 0.6731244 0.001 0.006 *
## 4 Metro vs Amp 16.937747 0.4241037 0.001 0.006 *
## 5 Metro vs Amp + Metro 24.788139 0.6077291 0.001 0.006 *
## 6 Amp vs Amp + Metro 27.745340 0.5691896 0.001 0.006 *
pairwise.adonis(ps1_otu_table.d20[,1:1363], ps1_otu_table.d20$Treatment)
## [1] "Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"
## pairs F.Model R2 p.value p.adjusted sig
## 1 Vehicle vs Metro 6.82161 0.2748254 0.001 0.006 *
## 2 Vehicle vs Amp 16.74809 0.4213559 0.001 0.006 *
## 3 Vehicle vs Amp + Metro 47.10367 0.7464490 0.001 0.006 *
## 4 Metro vs Amp 21.72066 0.4856963 0.001 0.006 *
## 5 Metro vs Amp + Metro 109.84593 0.8728604 0.001 0.006 *
## 6 Amp vs Amp + Metro 19.88933 0.4864185 0.001 0.006 *
# Amp ~ Day
ps1.amp.d0 <- subset_samples(ps1.amp, DaysTreatment == "D0")
ps1.amp.d3 <- subset_samples(ps1.amp, DaysTreatment == "D3")
ps1.amp.d7 <- subset_samples(ps1.amp, DaysTreatment == "D7")
ps1.amp.d13 <- subset_samples(ps1.amp, DaysTreatment == "D13")
ps1.amp.d16 <- subset_samples(ps1.amp, DaysTreatment == "D16")
ps1.amp.d20 <- subset_samples(ps1.amp, DaysTreatment == "D20")
doadonis(ps1.amp.d0, "SurvivalStatus") # *
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.05114 0.051141 1.0061 0.07183 0.418
## Residuals 13 0.66079 0.050830 0.92817
## Total 14 0.71193 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0040013 0.004001300 2.693166 999 0.117
## Residuals 13 0.0193144 0.001485723 NA NA NA
doadonis(ps1.amp.d3, "SurvivalStatus") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.15533 0.15533 0.83339 0.06024 0.647
## Residuals 13 2.42299 0.18638 0.93976
## Total 14 2.57832 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.000335951 0.000335951 0.01707187 999 0.917
## Residuals 13 0.255822190 0.019678630 NA NA NA
doadonis(ps1.amp.d7, "SurvivalStatus") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.13571 0.13572 1.1761 0.08297 0.25
## Residuals 13 1.50008 0.11539 0.91703
## Total 14 1.63579 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.007527834 0.007527834 2.554168 999 0.124
## Residuals 13 0.038314564 0.002947274 NA NA NA
doadonis(ps1.amp.d13, "SurvivalStatus") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.21715 0.21715 1.4577 0.10083 0.136
## Residuals 13 1.93658 0.14897 0.89917
## Total 14 2.15374 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.02576481 0.025764813 4.057142 999 0.075
## Residuals 13 0.08255628 0.006350483 NA NA NA
doadonis(ps1.amp.d16, "SurvivalStatus") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.20463 0.20463 1.0978 0.07787 0.328
## Residuals 13 2.42313 0.18639 0.92213
## Total 14 2.62776 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01037965 0.010379646 1.535046 999 0.213
## Residuals 13 0.08790315 0.006761781 NA NA NA
doadonis(ps1.amp.d20, "SurvivalStatus") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.1373 0.13734 0.58365 0.04297 0.887
## Residuals 13 3.0591 0.23531 0.95703
## Total 14 3.1964 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.001757289 0.001757289 0.1385605 999 0.695
## Residuals 13 0.164872085 0.012682468 NA NA NA
# Treatment ~ Day
ps1.d0 <- subset_samples(ps1, DaysTreatment == "D0")
ps1.d3 <- subset_samples(ps1, DaysTreatment == "D3")
ps1.d7 <- subset_samples(ps1, DaysTreatment == "D7")
ps1.d13 <- subset_samples(ps1, DaysTreatment == "D13")
ps1.d16 <- subset_samples(ps1, DaysTreatment == "D16")
ps1.d20 <- subset_samples(ps1, DaysTreatment == "D20")
doadonis(ps1.d0, "Treatment") # *
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 3 0.36198 0.120660 2.2241 0.13996 0.001 ***
## Residuals 41 2.22432 0.054252 0.86004
## Total 44 2.58630 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.006101085 0.002033695 1.53105 999 0.231
## Residuals 41 0.054460352 0.001328301 NA NA NA
doadonis(ps1.d3, "Treatment") # **
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 3 3.5719 1.19063 9.8536 0.41894 0.001 ***
## Residuals 41 4.9541 0.12083 0.58106
## Total 44 8.5260 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.3083364 0.102778794 11.34675 999 0.001
## Residuals 41 0.3713778 0.009057995 NA NA NA
doadonis(ps1.d7, "Treatment") # ***
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 3 3.1929 1.06429 11.477 0.45645 0.001 ***
## Residuals 41 3.8021 0.09273 0.54355
## Total 44 6.9949 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.3818544 0.127284804 33.01172 999 0.001
## Residuals 41 0.1580856 0.003855746 NA NA NA
doadonis(ps1.d13, "Treatment") # ***
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 3 4.1878 1.3959 11.498 0.45692 0.001 ***
## Residuals 41 4.9775 0.1214 0.54308
## Total 44 9.1653 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.3405377 0.113512572 14.31893 999 0.001
## Residuals 41 0.3250254 0.007927448 NA NA NA
doadonis(ps1.d16, "Treatment") # ***
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 3 5.3900 1.79667 14.644 0.52974 0.001 ***
## Residuals 39 4.7849 0.12269 0.47026
## Total 42 10.1749 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.3688410 0.122946985 17.11128 999 0.001
## Residuals 39 0.2802205 0.007185142 NA NA NA
doadonis(ps1.d20, "Treatment") # ***
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 3 5.2036 1.73452 15.837 0.54919 0.001 ***
## Residuals 39 4.2714 0.10952 0.45081
## Total 42 9.4750 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.8190444 0.273014799 41.90247 999 0.001
## Residuals 39 0.2541038 0.006515482 NA NA NA
sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2 pairwiseAdonis_0.0.1 cluster_2.0.6
## [4] mgcv_1.8-22 nlme_3.1-131 data.table_1.10.4-3
## [7] ggpubr_0.1.6 magrittr_1.5 microbiome_1.0.0
## [10] plotly_4.7.1 knitr_1.18 gridExtra_2.3
## [13] vegan_2.4-5 lattice_0.20-35 permute_0.9-4
## [16] RColorBrewer_1.1-2 phyloseq_1.23.1 plyr_1.8.4
## [19] reshape2_1.4.3 forcats_0.2.0 stringr_1.2.0
## [22] dplyr_0.7.4 purrr_0.2.4 readr_1.1.1
## [25] tidyr_0.7.2 tibble_1.4.1 ggplot2_2.2.1.9000
## [28] tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] lubridate_1.7.1 httr_1.3.1 rprojroot_1.3-2
## [4] tools_3.4.3 backports_1.1.2 R6_2.2.2
## [7] lazyeval_0.2.1 BiocGenerics_0.24.0 colorspace_1.3-2
## [10] ade4_1.7-10 withr_2.1.1 tidyselect_0.2.3
## [13] mnormt_1.5-5 compiler_3.4.3 cli_1.0.0
## [16] rvest_0.3.2 Biobase_2.38.0 xml2_1.1.1
## [19] labeling_0.3 scales_0.5.0.9000 psych_1.7.8
## [22] digest_0.6.14 foreign_0.8-69 rmarkdown_1.8
## [25] XVector_0.18.0 pkgconfig_2.0.1 htmltools_0.3.6
## [28] htmlwidgets_0.9 rlang_0.1.6 readxl_1.0.0
## [31] rstudioapi_0.7 shiny_1.0.5 bindr_0.1
## [34] jsonlite_1.5 crosstalk_1.0.0 biomformat_1.6.0
## [37] Matrix_1.2-12 Rcpp_0.12.14 munsell_0.4.3
## [40] S4Vectors_0.16.0 ape_5.0 stringi_1.1.6
## [43] yaml_2.1.16 MASS_7.3-48 zlibbioc_1.24.0
## [46] rhdf5_2.22.0 grid_3.4.3 parallel_3.4.3
## [49] crayon_1.3.4 Biostrings_2.46.0 haven_1.1.0
## [52] splines_3.4.3 multtest_2.34.0 hms_0.4.0
## [55] pillar_1.0.1 igraph_1.1.2 codetools_0.2-15
## [58] stats4_3.4.3 glue_1.2.0 evaluate_0.10.1
## [61] modelr_0.1.1 httpuv_1.3.5 foreach_1.4.4
## [64] cellranger_1.1.0 gtable_0.2.0 assertthat_0.2.0
## [67] mime_0.5 xtable_1.8-2 broom_0.4.3
## [70] survival_2.41-3 viridisLite_0.2.0 iterators_1.0.9
## [73] IRanges_2.12.0